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Abstract 

Similarity solutions play an important role in many fields of sci- 
ence. The recent book of Barenblatt j|] discusses many examples. 
Often, outstanding unresolved issues are whether a similarity solution 
is dynamically attractive, and if it is, to what particular solution does 
the system evolve. By recasting the dynamic problem in a form to 
which centre manifold theory may be applied, based upon a trans- 
formation by Wayne ju]], we may resolve these issues in many cases. 
For definiteness we illustrate the principles by discussing the appli- 
cation of centre manifold theory to a particular nonlinear diffusion 
problem arising in filtration. Theory constructs the similarity solu- 
tion, confirms its relevance, and determines the correct solution for 
any compact initial condition. The techniques and results we discuss 
are applicable to a wide range of similarity problems. 
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1 Introduction 

Consider the nonlinear diffusion problem with a step in the diffusivity dis- 
cussed by Barenblatt 0, §3.2] which in nondimensional form is 

a _ J 9 XX , 9 t >0 , , 

Ut -\ (l + e)0„, 9 t <0 ' ^ 

where 6(x,t) is the evolving concentration of some spatially distributed sub- 
stance. Such a problem, with its nonlinear step in the diffusivity, arises in 
theory of filtration of an elastic fluid in an elasto-plastic porous media (see the 
discussion in |2], §3.2.1]). It describes the diffusion in one spatial dimension 
x which is assumed here to be effectively of infinite extent. 

We write and analyse (|l|) as a perturbation of the basic linear diffusion 
problem, namely 

0t = xx + f(6,e), (2) 
where, since 9 t has the same sign as 9 XX , 

f _\0, 9 XX >0 ,g\ 
\ t9 xx , 9 XX < 

The term f(9, e) acts as a nonlinear perturbation to the basic diffusion of 

9 t = 9 XX (4) 

on an infinite domain. Of course e need not be small but we shall treat it so. 

We apply centre manifold theory to help understand and solve this prob- 
lem. But on the infinite spatial domain there is no clear cut centre eigenspace 
for (f|). However, following Wayne [l^, [| we transform the problem to one 
of seeking r) where 

r = logt, £=4> = 7f 0(r '^- (5) 

Then the dependence upon the scaled space variable £ causes the Gaussian 
spread from a point release, 
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to correspond to a fixed point of the dynamics for 0, namely 

a =e^ 2 / 4 . (7) 



Also, the algebraic decay in t from any compact release to the Gaussian (|6]) 
transforms to an exponentially quick decay in r to the fixed point (|?|). Centre 
manifold theory is applied in Section 2 to justify the self-similar Gaussian (||) 
as a valid approximation to the long-term dynamics of the non-constant 
diffusivity problem (|l|). Then the centre manifold analysis, as extended in 
Section 3, determines that the amplitude a of the decaying Gaussian evolves 
like 

a « a r €/V ^ (8) 

in accordance with the result reported by Barenblatt for e ^ 0. In addition 
to this confirmation of earlier results, centre manifold theory || immediately 
guarantees the attraction of the similarity solution. That is, this approach 
easily establishes the relevance of the similarity solution to the long-term 
dynamics of this nonlinear diffusion and we expect it to be able to analogously 
justify the relevance of similarity solutions for other problems. 

The amplitude of the spreading Gaussian not only decays in time, it also 
is a function of the initial distribution 6(x, 1) of the substance (note that the 
initial release is assumed to occur at t = 1 corresponding to the transformed 
time r = 0). Qualitatively, the long term behaviour is similar for all initially 
compact releases. However, the specific evolution of the model does depend 
on the specific initial conditions. In other words, we need to determine ao 
in (|3|). Naively we may expect that the total amount of substance in the 
model, given by a in (||), will be the same as that at the instant of release 
and so use _ 

roc 

a o — / 0{x, 1) dx . (9) 



However, this is only a leading order approximation and needs correction de- 
pending upon other details of the release distribution 6(x, 1). The corrections 
cannot be determined by scaling law arguments, but require a knowledge of 
the dynamics of approach to the similarity solution. Recently developed the- 
ory 0) is used in Section 4 to determine the proper choice of the initial 
conditions for the model amplitude a. 

For any given release of substance, the assumed origin of space-time may 
not be the best location for the origin of the similarity solution. In Section 5 
we show how the translational degrees of freedom in the coordinate system 
can be incorporated into the model for it to represent better the solution 
of the original diffusion problem. Numerical solutions reported in Section 6 
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confirm the effectiveness of the correct choice of as well as of time and 
space origins of the model. 

Finally we comment that the example discussed in detail here is just 
one of a wide class of nonlinear advection-reaction-diffusion problems. Cen- 
tre manifold theory may be successfully applied to many of these problems 
and not only create the similarity solution, but also justify its relevance as 
an attractive manifold, and determine the correct initial amplitude for the 
similarity solutions. One class of nonlinear reaction-diffusion problems was 
similarly analysed by Gene Wayne |10| . Some of the similarity solutions of 
the nonlinear advection diffusion problems discussed by Doyle and Englefeld 
II are also amenable to this centre manifold approach. 



2 Similarity solutions form a centre manifold 

Now investigate the centre manifold analysis in more detail. The transfor- 
mation (H) changes (Q) to 

r = £0 + /(0,e), (10) 

where the linear operator 
Adjoin the trivial equation 

e T = 0. (12) 

Then observe that for e = the Gaussian (0) describes a fixed point of fllQl)- 
(|T2| ) for all amplitudes a. Thus the centre manifold we construct will be 
global in a and local only in e. Now the linear operator L has a spectrum of 

X = -n/2, n = 0,1,2,.... (13) 

This is straightforwardly shown by looking for eigensolutions in the form 
e A„i--§ /4jy^£^ where H n are Hermite polynomials |]J. With two zero eigen- 
values, one from ( |i~3"l) and one trivially from (0), and the rest strictly neg- 
ative, centre manifold theory asserts there exists a two dimensional centre 
manifold for (|ll])-([l2|), A4 C , which is exponentially attractive to nearby tra- 
jectories. 

Thus by Theorem 2 in || p. 4], centre manifold theory immediately proves 
the attraction to the asymptotic similarity solution, albeit only for small 
enough e. (Contrast the ease of obtaining this result with Barenblatt's sta- 
bility analysis 0, §8.3.2].) In agreement with Barenblatt's equation (8.67), 
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from the spectrum (|T3|), we immediately deduce that the longest-lasting tran- 
sient in the approach to the similarity solution will be of relative magnitude 
approximately e^ T ^ 2 = 1/ \ft. 

We now approximate A4 C , parameterized by a and e, and the evolution 
thereon by 



= a(r) [-0o(O + eMO + (?M£) + O (e 3 )] , where 0> o 



e -? 2 /4 



2v^F ' 

(14) 

s.t. a = ag = a eg x + e 2 g 2 + O (e 3 ) (15) 

(■00 is normalised such that J^ipodC, = 1 and the overdot denotes d/dr). 
Substituting (|T4]) and (|15|) into (pi]) and equating all terms of O (e) we need 
to solve 

Cipi = ipo9i ~ D&ipo , (16) 

where for any s 

n I ' ^l" s - s l ,m 
D - = \&, ' (17) 

Here £o = is such that 0o^( — £o) = ^o^(^o) = 0. But C is singular as 
it has a zero eigenvalue; so we choose g\ to put the remaining terms in the 
range of £ — this is the solvability condition. In order to do this we take the 
inner product of equation ( [IB] ) with the solution z of the adjoint problem 

tfz = ZK-±Zzt = 0, (18) 
where the adjoint is obtained using the obvious inner product 

(u,v) = uvd£. (19) 



For a reason discussed later in the paper we normalise the adjoint eigenvec- 
tor such that (z,ipo) = 1. It is straightforward to check that the adjoint 
eigenvector satisfying this normalisation is z = 1. 

Finally, applying the solvability condition we find that 

g 1 = 2ifa{Z a ) = --L=. (20) 
V27re 

(As usual, we do not need to find ip\ to determine the leading order evolution.) 
The leading order centre manifold model a ~ —ea/\ / 2ne then has solution 



a = aoe 



-ir/V^Te = aot -a/2 ? where a= _ (21) 
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in agreement with Barenblatt [|2], pp 175-6]. The constant ao is determined 
by the initial conditions for the full original problem and will be determined 
in Section 4. 



The next-order correction matches earlier 
results 



Before proceeding to the next order approximation for the evolution on the 
centre manifold we need to find ipi. 

Since the operator £ is singular the solution is not unique and we are free 
to impose one additional condition on the solution to fix it. It is convenient 
to require that 

poo 

= . (22) 



Physically this implies that the total amount of the diffused substance is 
given completely by the leading order approximation of the solution, and 
as ff^ ipo = 1, the total amount is simply a. Under this condition the 
continuous, up to the second derivative, solution to fll6"D becomes 



■0i 



1 C 3 + ^ (erf(f) - l) erf(f 
/«erf(f)e^ 2 / 4 ^ 



8V7T ~"~ 2\ 2« 



(erf(f)-erf(-) 
where H denotes the Heaviside function and 



(23) 



C 3 



27rV2e 



(24) 



The integrals entering the definition of C3 are: 



h 



6) 



e 4 erf - erf — 



/ 



6> 



erf U - 1 



erf f — 



di 



-0.1076980691 



0.2262196880* 



-0.1358229603i 



e 4 







e 4 erf 







dyd£ » 0.6931471806i . 



(25) 

(26) 
(27) 



As expected the first order correction, ^1, is an even function of £, see Fig- 
ure [TJ. 
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0.2 - 



0.1- 








-6 -4 -2 2 4 6 

Figure 1: Solutions ?/>o(0 (solid line) showing the Gaussian shape of the 
basic similarity solution, and (dashed line) showing that the Gaussian 

is flattened and broadened by the nonlinear diffusion. 
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Let (£) = 0. Then £ = £o + e£i + (e 2 ) where, as is deduced from (|T 
and (H), 

A = - « 0.5665706981. (28) 

V><*«(£o) 

Collecting terms of (9 (e 2 ) we obtain 

£ip2 = i>i9i + 1P092 - (-Dft+efi - ^o) - D&ifo. . (29) 
Similarly to the previous section, the application of the solvability condition, 



upon making use of fl22|) , leads to 



g 2 = 2 (^(£ + e£i) + Vo^Co + efi) - ^Mfo)) 

= 2^(£ o ) + 0(e) (30) 
« 0.06354624322 + O (e) , 

where the even symmetry of ?/>o and ?/>i is taken into account. The numerical 



results given in ( p8|) and ( j30|) coincide with the ones reported by Cole and 
Wagner in their paper j|, p. 167] though our values are given with more 
significant digits. Consequently, the next order centre manifold model is 

a « a(e#i + e 2 g 2 ) (31) 

with solution 

a = a r a ' /2 , where a' = 2e ( -== - eg 2 ) . (32) 

\V27re / 

4 The correct initial condition ensures fidelity 
of the model 

The correct projection of initial conditions onto a centre manifold, first de- 
veloped in and recently refined in M , should approximately determine the 
"functional of the initial conditions" mentioned by Barenblatt near the top of 
p. 202 [Q, but not previously found. Here we follow the procedure outlined in 
|| to give the proper initial conditions a for the centre manifold model (|32|) 
when the initial conditions for the original problem are given by 9 = 9q(x) 
at t — 1 corresponding to r = 0. We expect that a\ T=0 = J^ oo 9 dx, but 
this is only a first approximation. The more careful analysis corrects this 
approximation. 

As used in previous sections, the special form of (|l0f) implies that its 
solution is to be found in the separable form 



</>(t, £; e) = a(r)^(C; e), where a = a(r)g(e). (33) 
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Then "vectors" locally tangent to the centre manifold are found to be ei = 
(adip/de,l) and e 2 = (ip,0). According to [§] we need to find "vectors" z 1 
and Z2 satisfying 



T>Zj - ^2 (Dzj, e k ) z fc = , j = l,2 



(34) 



k=l 



and normalisation (zj, e&) = 5jk where the dual operator T> is defined as 



the adjoint 



V 



O + eD\ 








and 



d 



D\ = D- i + 2 (6(£ + - 6(£ - 0) ^ + 6'(Z + " " 



(35) 



(36) 



(37) 



in which 5 and 5' denote the Dirac delta function and its derivative, respec- 
tively. The normalisation conditions give that 



(1) 
2 

(2) 



0. 
1. 



roo / _i .. 

J-OO I ' 1 07 ' 



/- 



oo V 2 Se 



(2) 
1 

,(2) 



arf (0 

de = i 



+ ^ = 



(38) 
(39) 



We look for the solution of (|34]) satisfying Pzi = 0, i.e. 













ri 



(40) 



Hence we immediately deduce that = 0. Consequently, the second of 
normalisation conditions (|3~9"D is transformed to J^r^dC, = 1. Then from 
the projection of initial conditions 



— (r?, 0o - a| r=0 ^) + (eo - e) (r? \ l) 







a| T= o 



(41) 



and we deduce that e = e . This result, that the parameter e remains un- 
changed between the model and the original problem, is expected at the 
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outset, but we have just demonstrated how it is obtained in the context of 
the developed theory for the projection of initial conditions. 

Thus the proper initial condition for the amplitude a\ T= o is given by 

(4 1) ,^ -a| r= oV) = 0, (42) 

or, equivalently, since the problem is linear in amplitude a and the normali- 
sation conditions (|39"D are used, by 



a\ T =o 



(4 1] ,9 ). (43) 



Thus the problem of finding the proper initial condition is reduced to solving 
for which satisfies the following equation deduced from (|34|) 



(£t + 6 Dt) 4 1] = ((& + eDt) rf } , rf 3 . (44) 



Performing integration by parts in the right-hand side of (f44|) and using the 
normalisation (|39|) we obtain 

(£ f + eD\) - gr? = 0, (r?\ ^) = 1- (45) 

We solve (|4"5| ) assuming r% = Po(0 + e Pi(0 + O (e 2 ) and recollecting that 
g « -e/^/^ne + O (e 2 ) and ^ = ifo + e^i + O (e 2 ). At (e°) we obtain 

tfpo = 0, (po,^o) = l (46) 

with solution p — z — 1. Thus at leading order a| T=0 = /f^ #o(£)^£- 
At C (e 1 ) we obtain 

^ + + 6) - S'(£ - Co) + -i= = 0, (ft, Vo> = 0. (47) 

V27re 

The solution, presented in Figure 0, has the following algebraic form 

Pl (0 = c 4 + (l + z^/Jerf (^)) (JT(£ - - # (£ + &)) 

-vti"oerf(|)e^ 2 /4 dy (48) 

+*vi ( x + erf (i) - F ^ + - H (£ - &)) erf ( 2 ) > 



where 



c 4 = 7=^3 - v^(/ 2 + h)) + (l + »7Eerf (^)) erf (A 
« 0.0589390531. 



(49) 
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Finally we then have that the proper initial condition for the centre man- 
ifold model (^l|) is given by 

a\ T=0 = /_"(! + epi(0 + O (e 2 )) 8 (f) d£ (50) 

Note that p 1 ~ [2/ (7re)] 1//2 log(|,f|) as |£| — > oo and, consequently, the inte- 
gral ( f43|) converges only for a sufficiently compact initial distribution 6*0. This 
emphasises that the projection of the initial conditions is local in its nature 
and it is applicable only if the initial conditions for the original problem are, 
in some sense, close to the centre manifold. 



5 Choose an optimal origin in time and space 

It follows from the transformation of space and time variables (|5|) that the 
diffusion from a localised initial release of arbitrary form occurring in the 
original problem at t = 1 is modelled by the evolution from the initial state of 
a point release, a delta function, at x = t = 0. On the other hand the original 
partial differential equation (p]) is invariant with respect to translations in 
time and space. Thus there is freedom to choose the time and space origins 
for the model to suit best the actual distribution of the initial 6. To account 
for these inherit degrees of freedom in the original problem we generalise the 
coordinate transformation (|]) to 

r = log(t + t ), ^ = 4^, e=*}££L, (51) 

where to > 0. Now the localised release 9q(x) occurring in the original prob- 
lem at time t — (not at t — 1 as assumed in the previous sections) is 
modelled by some Gaussian centred at x$ rather than by the delta function 
at x = 0. The width of the model Gaussian at the moment of the actual 
release t = is determined by to which also determines the location of the 
virtual origin in time for the model. Generalisation (EWj does not affect the 
analysis of the previous sections. In particular, the model dynamics (|31|) 
is unchanged because the general long-term dynamics are independent of 
the space-time origin. However, the generalisation provides a two-parameter 
family of model solutions to the original problem ([]]) rather than just the 
unique model described earlier. Thus here the general projection of initial 
condition (EDI) becomes 



/■oo 

a 



6 Q (x)dx. (52) 
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One is free to choose parameters Xq and to entering (^) in such a way that 
the model possess certain additional properties. For instance, we choose t 
such that the contribution of the e-dependent terms in fl5"2|) is zero — this 
choice should ensure that the model a most closely matches the solution 9 
for the original problem in the short-term as well as the long-term evolution. 
In essence this is equivalent to considering all the centre manifolds (in a 
and e) parameterized by to and xq, and choosing that centre manifold whose 
isochrons are linearly "vertical" and hence make the definition of a match 
the projection. It is always possible to make this choice since physical initial 
distributions 9q are non-negative functions while the mean of p\ is zero. Thus 
require 

pi ( J 9 (x) dx = (53) 



oo 



/to / 

which we view as implicitly defining t as a function of Xq. 

The value of xq is then fixed to minimise to- We feel this is desirable 
since it minimises the spread of the model's Gaussian at the initial instant 
of release and so maximises the information content of the model. (It is also 
the only distinguished x .) Differentiating fl53|) with respect to x we obtain 



dl 1 f°° , / x-x \ 



dXQ y/i J-oo \ -\/t() 



X — Xn dtn 

1 + 



2t dx 



dx = 0, (54) 



where prime denotes differentiation with respect to the argument. At the 




point of extremum dto/dxQ = and the second term in the brackets in (54) 
vanishes. Thus we solve 

x) dx = . (55) 

in conjunction with ([53]) to define xq and to- As an aside it follows from 
the above discussion that such chosen xq and to guarantee that / = is a 
minimum contribution to the e-correction of initial conditions for the model. 
If 8q is symmetric, say about x = q, then, owing to the even symmetry of 
Pi, the choice of Xo = q guarantees that ( |55D is satisfied. Thus for symmetric 
6*o the best choice for the centre of the Gaussian spread of the model is the 
point of symmetry. 

Finally, the initial amplitude is then given by 



a 



t^ 72 r 9 Q (x)dx (56) 



and the model solution written in the original variables becomes 

(57) 



9 



(t + t )( 1+a ')/ 2 



6 NUMERICAL RESULTS DEMONSTRATE THE ACCURACY 



14 



where £ an d %o satisfy (^3|) and (|55|). 

6 Numerical results demonstrate the accu- 
racy of the model 

We illustrate the correctness of the derived initial conditions by comparing 
the model predictions with the direct numerical integration of equation ([[J). 
Let the initial distribution of substance for the original problem at t = be 
in the form of the Gaussian 

9 Q = W — exp(-10x 2 ) . (58) 

V 7T 

Numerical integration of ([!]) with initial distribution (|58 ) was performed 



using IMSL routine dmolch with the accuracy of 10~ 8 . Since the long 
term behaviour of the numerical solution was found to depend on the size 
of the computational domain, the preliminary test of the numerical solution 
was performed for e = for which the analytic solution comes from It 
was found that the non-reflecting boundary conditions 6 X (L)/8(L) = x/(2t) 
imposed at L = 22.5 eliminated such an influence for the time interval con- 
sidered. 

The resulting time evolution of the direct numerical solution for e = 0.1 
at x = is shown by a solid line in Figure Because of the symmetry of 
initial distribution fl58|) with respect to the line x = 0, (|^) gives the value 
Xq = for model ((57|). Numerical evaluation shows that condition (|53| ) is 
satisfied for 8 given by flSlf ) for t rs 0.0250. As seen from Figure |3](a) the 
model dynamics shown by star symbols approaches the numerical curve very 
quickly. In Figure |3](b) we compare the numerical and the proper model (|57|) 
solutions with the earlier proposed model B. [1 



(59) 



which uses naive initial condition (g) — shown by diamond symbols. While 
the present model and numerical solution are virtually indistinguishable in 
their long term evolution, model ( |5T?D based solely on scaling arguments is 
able to predict just a slope. The actual values of the distribution maximum 
it provides lies apart from the numerical curve for all time. Thus the correct 
initial conditions for the model are essential to avoid a permanent finite phase 
difference between the model and the actual full solutions. 
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-4 -2 2 1.0 1.5 2.0 2.5 3.0 

log(t) log(t) 

Figure 3: Numerical (solid line) solutions of equation (P evaluated at x — 
for e = 0.1 compared with the model (0) that uses the correct initial 
conditions (stars) and the previous model (|59D (diamonds). 
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7 Conclusions 

We have demonstrated that the centre manifold theory provides a straightfor- 
ward and rigorous way of deriving not only the functional form of similarity 
solutions of nonlinear diffusion, but also the appropriate initial conditions for 
the model in terms of the initial distributions of the substance. This cannot 
be done using other modelling approaches such as, for example, scaling laws 
or the method of multiple scales. The correct provision of initial conditions 
also enables us to determine an optimal location for the virtual space-time 
origin for the model. The present technique may be successfully used for 
modelling a wide class of nonlinear filtration/diffusion problems. 
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